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A diffusion model of the time evolution of loss rates caused by a step in collimator position 
is presented. It builds upon the model of Ref. [1] and its assumptions: (1) constant diffusion 
rate within the range of the step and (2) linear halo tails. These hypotheses allow one to 
obtain analytical expressions for the solutions of the diffusion equation and for the corre- 
sponding loss rates vs. time. The present model addresses some of the limitiations of the 
previous model and expands it in the following ways: (a) losses before, during, and after the 
step arc predicted; (b) different steady-state rates before and after are explained; (c) deter- 
mination of the model parameters (diffusion coefficient, tail population, detector calibration, 
and background rate) is more robust and precise. These calculations are the basis for the 
measurement of transverse beam diffusion rates as a function of particle amplitude with 
collimator scans. The results of these measurements in the Tevatron will be presented in a 
separate report. 
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I. 



INTRODUCTION 



Phenomena related to the stochastic transverse beam dynamics in circular accelerators can 
be described in terms of particle diffusion [2-6]. It was demonstrated that these effects can be 
observed with collimator scans [1]. Usually, collimator jaws are the devices that are closest to the 
beam and they define the machine aperture. If they are moved towards the beam center in small 
steps, typical spikes in the local shower rate are observed, which approach a new steady-state level 
with a characteristic relaxation time. When collimators are retracted, on the other hand, a dip in 
losses is observed, which also tends to a new equilibrium level (Figure 1). A detailed description 
of the Tevatron collimation system can be found in Ref. [7]. 

These phenomena have been used to estimate the diffusion rate in the beam halo in the SPS at 
CERN [8], in HERA at DESY [1], and in RHIC at BNL [9]. Similar measurements were carried 
out at the Tevatron in 2011. Besides the interest in characterizing the beam dynamics of colliding 
beams, these measurements were motivated by the study of the effects of the novel hollow electron 
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FIG. 1. Example of loss rate data taken during a collimator scan in Tevatron Store 8749 (20 May 2011): 
local loss rate (top, device T:LF480); collimator position (bottom, device T:F48VCP). The collimator steps 
take about 0.2 s. In this example the collimator position was recorded only every 15 s. 



4 



beam collimator [10]. 

Here we present a more complete model of beam evolution under diffusion. It will serve as the 
basis for interpreting of Tevatron data. Previous models are extended to explain the behavior of 
losses before, during, and after the collimator step. This allows one to extract the diffusion rate in 
a more robust way, by taking into account not only the relaxation time, but also the steady-state 
loss rates before and after the step and the peak or dip value. The analysis of Tevatron data will 
be presented in a separate report. This model can also be applied to the dynamics of beams in the 
LHC. 



II. MODEL 



Following Ref. [1], we consider the evolution in time t of a beam of particles with phase-space 
density f(J,t) described by the diffusion equation: 

d t f = dj(Ddjf) (1) 

where J is the Hamiltonian action and D(J) the diffusion coefficient. The particle flux at a given 
location J = J' is (f> = —D • [djf] J=J i. 

During a collimator step, the action J c = x\j ' j3 c , corresponding to the collimator position x c at 
a ring location where the amplitude function is /3 C , changes from its initial value J c j to its final 
value J c f during a time At. The step in action is A J = J c f — Jd. In the Tevatron, typical steps 
are 50 fim. in 0.2 s, and the amplitude function is tens of meters. The behavior of J c (t) can be 
modeled, for instance, by a linear function connecting J c j with J c f. 

Jd t<0 
Jc{t) = I J ci + (Jcf - Jd) ■ t/At < t < At (2) 
Jcf At < t 

It is assumed that the collimator steps are small enough so that the diffusion coefficient can 
be treated as a constant in that region. This hypothesis is justified by the fact that the fractional 
change in action is of the order of AJ C /J C ~ (2) (25 /xm)/(2 mm) = 2.5%. Because the diffusion 
coefficient is a strong function of action (D ~ J 4 ), this translates into a variation of 10% in the 
diffusion rate, an acceptable systematic in a quantity that varies by orders of magnitude. If D is 
constant, the diffusion equation becomes 



dtf = Ddjjf. 



(3) 



5 



With these definitions, the particle loss rate at the collimator is 



L = -D.[djf) J=Jc . 



(4) 



Particle showers caused by the loss of beam are measured with scintillator counters placed close to 
the collimator jaw. The observed shower rate is parameterized as follows 



S = kL + B, 



(5) 



where A; is a normalization constant including detector acceptance and efficiency and B is a back- 
ground term which includes, for instance, the effect of residual activation. Both k and B are 
assumed to be independent of collimator position and time during the scan. 

III. BOUNDARY CONDITIONS 



The collimator is treated as a perfect absorber, so that the boundary condition for the phase- 
space density becomes 



/(J, t) = for J > J c 



(6) 



c 

u 
c 

a 

c 

o 

I 




0.00 0.01 0.02 0.03 0.04 0.05 



Action [urn] 



FIG. 2. Illustration of boundary conditions in the case of an outward collimator step: initial distribu- 
tion fo{J), intermediate asymptotic distributions / a (J, t), final distribution foo{J)- The vertical lines repre- 
sent the positions of the collimator vs. time. Parameters in this example are J c .; = 0.04 fim, J c f = 0.05 (im, 
Ai = 1 pr 2 , A f = 0.8 /mi~ 2 , At = 1 s. 



We assume that cancellation of the particle flux at J = is automatically satisfied by D(0) ~ 0, so 
no boundary condition is imposed there. This greatly simplifies the form of Green's function (see 
below) . 

As initial conditions and asymptotic behavior for the phase-space density we use linear functions 
of action: 



where A{ and Af are constants. This is the essential hypothesis that allows one to obtain analytical 
solutions for the time evolution of the distribution function. It is justified by considering these 
expressions as the first term in the Taylor expansion of the beam tails. A linear behavior of the 
asymptotic solution also gives a constant steady-state flux. 

In this respect, the present model differs from that of Ref. [1]. We allow the slopes of the initial 
and final distributions to be different. This is necessary to explain the difference in the steady-state 
loss rates L before and after the collimator step. If one assumes the same diffusion coefficient and 
the same slope before and after, then one can only predict the same steady-state rate. Including 
the measured steady-state loss rates before and after the collimator step helps to disentangle the 
effects of population and diffusion, and to give a physical meaning to the model parameters. 

Because of its linearity, a solution of the diffusion equation can be found using the method of 
Green's functions: 



where G(J,J',t) is Green's function for the given problem, /o (Eq. 7) is the initial distribution, 
and f a is an asymptotic solution. We are looking for a model of losses not only before and after 
the collimator step, but also as the collimator is moving. In this respect, too, this model extends 
that of Ref. [1]. For this reason, we use f a = foo (Eq. 8) for J c = J c f (i.e., t > At) and 




(7) 




(8) 




(9) 




(10) 
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as the collimator moves (0 < t < At). The parameter A c (t) is chosen to vary linearly between Ai 
and Af, 



A c (t) 



A, 



t < 



Ai + {A c - Ai) • t/At < t < At 
A f At<t 



(11) 



so that the asymptotic solution transitions smoothly from Eq. 7 to Eq. 8. The initial and asymptotic 
solutions are illustrated in Figure 2. 

The basic kernel for the diffusion equation is 



K(J,j',t) 
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(12) 



with a = y/2Dt. To satisfy the boudary condition at the collimator (Eq. 6), an antisymmetric 
Green's function can be used: 



G(J,J',t) = {exp 



1 f(J c -J')-(Jc-J)\ 



a 



exp 



1 gj c - j') + (j c - j) 

2 V a 



2no 



(13) 



so that 



G(J = J c , J', t) = G(J, J' = J c , t) = 0. 



(14) 



The requirement that the solution be zero beyond the collimator position, 



G(J, J',t) = if J c < J or J c < J', 



(15) 



is automatically satisfied by limiting the integration region between and J c (Eq. 9). Imposing 
additional boundary conditions at J = would require G to be an infinite series. The analytical 
approximation used here does not constrain the phase-space density or its gradient at the origin, 
but it turns out a posteriori that f(0,t) does not vary significantly if /o(0) — /oo(0), which is what 
one would expect for a collimator step affecting the beam halo and not the beam core. Green's 
function also satisfies the general symmetry property G(J,J',t) = G(J' , J,t). We also note its 
asymptotic behavior in time: G(J, J',0) = 5(J — J') and G(J,J', oo) = 0, which justifies the 
physical interpretation of /o and f a in Eq. 9 as initial and asymtotic solutions. 
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IV. SOLUTIONS 

By setting up the diffusion model in the way described above, solutions can be expressed analyti- 
cally through Eq. 9. It is convenient to treat the cases of inward ( J c f < J c j) and outward (J c j < J c f) 
movement separately. In the inward case, the integrand is fo — f a = Ai(J c i — J) — A C (J C — J). In 
the outward case, it is convenient to divide the integral into two parts: 
rJc fJd rJc 

/ = / + / • (16) 
Jo Jo Jj ci 

This is done because /o is null beyond the initial collimator position: fo — fa = —fa ( see a ls° 
Figure 2). To express the primitive of the Gaussian function, we use the cumulative Gaussian 
distribution function P(x), defined in Appendix A. (Another possible choice is the so-called error 
function.) Integration yields the solutions of the diffusion equation in the two cases, fi(J, t) (inward 
step) and fo{J-,t) (outward), subject to the boundary conditions specified above: 

fi(J,t)= - Ai(J ci + J-2J c ) (17) 
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FIG. 3. Evolution of distribution function during collimator step: fi(J,t) (inward, left) and fo(J,t) (out- 
ward, right). The vertical lines represent the positions of the collimator vs. time. Collimator action varies 
between J C i = 0.05 /xm and J c f = 0.04 (im in the inward case (viceversa in the outward case) in a time 
At = 1 s. The slopes of the tails arc Aj = 0.8 /im" 2 and Af = 1 /im~ 2 in the inward case (viceversa 
outwards). The diffusion coefficient is D = 10~ 5 fim 2 /s. 
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Some examples of the evolution of the phase-space density described by these functions are shown 
in Figure 3. A few representative snapshots in time are chosen: during collimator movement; a 
short time after the step, with a time scale determined by t s = \J c i — Jcf\ 2 /D = 10 s; and a long 
time after the step, with a characteristic time ti = [min (J c i, J c f)] 2 /D = 160 s. 



V. TIME EVOLUTION OF LOSSES 



Local losses are proportional to the gradient of the distribution function at the collimator (Eq. 4) . 
The partial derivatives of the phase-space density with respect to action are the following: 



djfi(J,t) 



+ 



Ai + (Ai - A 
1 



2Ai(J c - J ci )exp 



1 (J c -J^ 
2 



a 



+ 









1 


exp 




+ exp 










~2 



a 



(19) 
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The value of the gradient at the collimator is therefore 

-Jc 



djfi(J c ,t) 



Ai + 2(A, - A C )P 



a 



(21) 
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These are the functions that are used to model the measured shower rates (Eqs. 4 and 5). As 
expected, both functions tend to — Ai for t — > and to —Aj as t — > oo. For the t — > limit to 
hold for the outward solution, it is necessary that J c — > J c i faster than \/t, which is satisfied by 
the linear approximation adopted here; otherwise, the slope will tend to zero. 

These functions explain the data very well. In the transient region, after the collimator has 
reached its final position, they agree with those calculated in Ref. [1]. Their main feature is a 
decay proportional to the square root of time (through the parameter cr), as is typical of diffusion 
processes. A few examples are plotted in Figure 4. 
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FIG. 4. Calculated evolution of loss rates L(t) during a collimator step according to Eqs. 4, 21 and 22: 
inward (left) and outward (right). Collimator action varies between J c ; = 0.05 /im and J c / = 0.04 fim in the 
inward case (viceversa in the outward case) in a time At = 1 s (see also Figure 3). The effect of 3 different 
values of the diffusion coefficient D is shown. The slopes of the tails are scaled so that the initial and final 
steady-state loss rates are the same in all cases: Ai = 1/D. Af = AiJ C i/J c f. 
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VI. COMMENTS ON PARAMETER ESTIMATION 

Having a model that describes the data before, during, and after the collimator step has several 
advantages. The products kDAi + B and kDAf + B are determined by the steady-state loss rates. 
If a data set includes measurements of several steps at different amplitudes, the parameters k 
and B, which are independent of J and t, can be determined separately. Therefore, steady-state 
rates constrain the products DAi and DAf at each step. 

The value of the diffusion coefficient D is constrained both by the peak (or dip) value relative 
to the steady-state rate and by the duration of the transient through the parameter a. In fact, 
the peak (or dip) value of the loss rate is achieved when the collimator reaches its final position 
(J c = J c f, t = At). At this point, neglecting the background, the loss rate is 

|AJ| 



S(At) ~ kDAi 



1 ± 



(23) 



y/irDAt. 

whereas S(0) — kDAi. It follows that an estimate of the diffusion rate is 

D ~ ^ (24) 

vrAi [S(At)/S(0) - l] 2 V ' 

On the other hand, losses relax with a typical time constant that depends on the diffusion rate and 
on the magnitude of the step. One may define the characteristic time t n so that 

l*£L * J= = 0.56, (25) 

meaning that the magnitude of the transient at time t n is about half of that of the steady-state 
rate. Therefore, we have an independent estimate of D: 

D^. (26) 

If only the t > At data is considered, as is done in Ref. [1], this is the only available information 
on D. In addition, in this case, the diffusion coefficient is highly correlated with the steady-state 
parameters. 

These rough estimates of the model parameters can be used as initial guesses in a least-squares 
fit of the data. 



Appendix A: Useful formulas 

To express the solutions of the diffusion equation, we use the cumulative Gaussian distribution 
function P(x), defined as follows: 



P(x) = -= I exp (-z z /2) dz. (Al) 
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Therefore, P(— oo) = 0, -P(O) = 1/2, and P(oo) = 1. By definition, the integral of a Gaussian 
function can be expressed as follows: 
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Appendix B: Scripts 

Numerical calculations, data analysis, and graphics were done with the open-source, multi- 
platform statistical package R version 2.12.0 (2010-10-15) [11]. This documentation was produced 
by integrating F>TpjX with R using the Sweave package. The source code can be found in the file 
dmcs.tar.gz. Below is the R part of the code. 

#line 71 "dmcs.Rnw" 

Stangle ('dmcs . Rnw' , annotate=FALSE) # make dmcs.R for later inclusion 



# collimator position vs. time, linear 

Jc <- function(t, Jci=0.04, Jcf=0.05, dt=l){ 

Jc <- 0*t + Jci # initial value 

subset <- dt<=t 

Jc [subset] <- Jcf # final value 
subset <- 0<t & t<dt 

Jc [subset] <- Jci + ( Jcf -Jci) *t [subset] /dt # collimator moving 
return( Jc) 

} 



# slope of distribution vs. time, linear 
Ac <- function(t, Ai=l, Af=0.9, dt=l){ 

Ac <- 0*t + Ai # initial value 

subset <- dt<=t 

Ac [subset] <- Af # final value 
subset <- 0<t & t<dt 

Ac [subset] <- Ai + (Af-Ai) *t [subset] /dt # collimator moving 
return(Ac) 
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# initial, final, and asymptotic distributions, linear 
flin <- functional, A=l, J0=0.04){ 

f <- 0*J # default outside range 

JO <- abs(JO) # define edge as positive 

subset <- abs(J)<J0 # triangular with cusp at J=0 

f [subset] <- A* (JO - abs (J [subset] ) ) # linear 

return(f) )• 

# evolution of the distribution function (der=0) 

# or its derivative (der=l) , analytical formulas 

fs <- function(J, t, Ai=l, Af=0.9, Jci=0.04, Jcf=0.05, D=le-5, dt=l, der=0) { 
if(der!=0 ft der!=l) stop ("Derivative ('der') must be or 1.") 
fs <- # default outside collimators 

if(Jci>Jcf) IN<-TRUE else IN<-FALSE # inward or outward step? 
if(t<0){ # initial distribution 

if(der==0) fs <- flin(J, A=Ai , J0=Jci) 

if(der==l) fs < Ai 

} else{ # evolution 

Jco <- Jc(t, Jci=Jci, Jcf=Jcf, dt=dt) # collimator position 

Aco <- Ac(t, Ai=Ai, Af=Af, dt=dt) # slope of asymptotic solution 

if(0<=J & J<=Jco){ # inside collimators 
s <- sqrt(2*D*t) # Gaussian sigma 

# terms common to both inward and outward step 

if(der==0) fs.base <- -pnorm(-J/s) * (Ai*(Jci-J) - Aco*(Jco-J)) + 
pnorm((J-2*Jco)/s) * (Ai* ( Jci+J-2*Jco) - Aco*(J-Jco)) + 
s/sqrt(2*pi) * ( 

exp(-0.5*(J/s) ~2) * (Aco-Ai) + 
(-l)*exp(-0.5*((2*Jco-J)/s)~2) * (Aco-Ai) ) 
if(der==l) fs.base <- (Ai-Aco) * (pnorm(-J/s) + pnorm( ( J-2*Jco) /s) ) + 
(Ai*Jci-Aco*Jco)/sqrt(2*pi)/s * 

(exp(-0.5*(J/s)~2)+exp(-0.5*((J-2*Jco)/s)~2)) 

# add specific terms 
if (IN) { # inward step 

if(der==0) fs <- fs.base + 

(-l)*Ai*(Jci+J-2*Jco) + 

2 * pnorm( ( Jco-J) /s) * Ai * (Jci-Jco) 
if(der==l) fs <- fs.base + 

(-l)*Ai + 2*Ai*(Jco-Jci)*exp(-0.5*((Jco-J)/s)~2)/sqrt(2*pi)/s 
}■ else { # outward step 
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if(der==0) fs <- fs.base + 

pnorm( (Jci-J) /s) * Ai * (Jci-J) + 

(-1) *pnorm( (Jci+J-2*Jco)/s) * Ai * ( Jci+J-2*Jco) + 
s/sqrt(2*pi) * ( 

exp(-0.5*((Jci-J)/s)~2) * Ai + 
(-l)*exp(-0.5*((Jci+J-2*Jco)/s)~2) * Ai ) 
if(der==l) fs <- fs.base + 

(-l)*Ai*(pnorm( (Jci-J) /s)+pnorm( (Jci+J-2*Jco)/s) ) 

} 

} 

} 

return(f s) 

} 

fv <- Vectorize(f s, vectorize . args=c ("J" , "t") ) # vectorize the function 

# loss rate vs. time 

# particle flux at the collimator L = -D * (df/dJ)|J=Jc 

# with conversion constant k and background B 

loss. rate <- function(t, Ai=l, Af=0.9, Jci=0.04, Jcf=0.05, D=le-5, dt=l, 
k=l, B=0) { 
L <- -Ai + 0*t # default before step 
subset <- t>=0 # evolution 
s <- sqrt (2*D*t [subset] ) # Gaussian sigmas 

Jco <- Jc (t [subset] , Jci=Jci, Jcf =Jcf , dt=dt) # collimator positions 
Aco <- Ac (t [subset] , Ai=Ai , Af =Af , dt=dt) # slopes 
if(Jci>Jcf) IN<-TRUE else IN<-FALSE # inward or outward? 
# common terms 

L [subset] <- 2* (Ai-Aco) *pnorm(-Jco/s) + 

2/sqrt(2*pi)/s * (Ai*Jci-Aco*Jco) * exp(-0.5*(Jco/s)~2) 
if (IN) { 

L [subset] <- L [subset] - Ai - 2/sqrt (2*pi) /s * Ai * (Jci-Jco) 
} else { 

L [subset] <- L [subset] - 2 * Ai * pnorm( (Jci-Jco) /s) } 
return(-k*D*L+B) 

} 

# sample data 

# read collimator position and loss rate data 
ReadACNETData <- f unction(fn="CFBI_NG • txt .gz") { 
cat ("Reading file ", fn, "...\n") 
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dn <- gsub(" .txt .gz" , "", fn) # name of data frame 

tO <- scan(fn, what =" char act er " , skip=l, nlines=l) # get initial time 
tO <- tO [6] 

h <- as .numeric (substr (tO , 1 ,2) ) 

m <- as .numeric (substr (tO ,4, 5) ) 

s <- as .numeric (substr (tO ,7 ,8) ) 

assignO'TO", h+m/60+s/3600 , envir= . GlobalEnv) 

cat ("Start time: ", TO, 11 h\n") 

cnam <- scan(fn, what="character" , skip=2, nlines=l, sep="\t") # column names 
cnam <- substr (cnam, 3, 12) # strip logger names 

cnam <- gsub(" [" [:alnum:]] " , "", cnam) # remove non alphanumeric characters 
if (length(cnam) <72) cnam <- gsub("0", "", cnam) # remove Os from scalars 
cnam [seq(l , length(cnam) ,by=2)] <- 

paste("t", cnam[seq(2,length(cnam) ,by=2)] , sep="") # times 
assign(dn, read. table (fn, header=FALSE , skip=3, colClasses="numeric" , 
na. strings=c("NA" , "999999") , col.names=cnam) , 
envir= .GlobalEnv) # fill data frame 

} 

ReadACNETData('Store8749_CP . txt . gz') 
ReadACNETData ( 'Store8749_TLF480_15Hz . txt . gz' ) 

tc <- Store8749_CP$tF48VCP/3600 + TO # time stamps of coll. positions [h] 

cpmil <- Store8749_CP$F48VCP # collimator position [mils] 

cpO <- 366 # collimator offset (arbitrary) 

cp <- (cpmil - cpO) * 25.4 # collimator position [urn] 

tTl <- Store8749_TLF480_15Hz$tLF48/3600 + TO # time stamps of losses [h] 
Tl <- Store8749_TLF480_15Hz$LF48 # loss monitor [Hz] 

# plot setup 

options (SweaveHooks=list (f ig=f unctionO { # common plot options 
pt<-10; fo<-"Times" 
ps . options (pointsize=pt , f amily=f o) 
pdf . options (pointsize=pt , f amily=f o) 
par(oma=rep(0,4) , mar=c(4. 2 ,4. 2 ,0 . 2 ,0 . 2) ) })) 

#line 217 "dmcs.Rnw" 

ti <- 9.105 # initial time [h] 

tf <- 9.294 # final time 
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plot. unit <- 60 

xl <- c(0, tf-ti)*plot.unit 

par(mar=c(0, 3.2, 0, 0), oma=c(3.2, 0, 0.2, 0.2), mgp=c (2 , 1 , 0) , tcl=-0.25) 
layout (matrixd : 2, nrow=2, ncol=l)) 

x <- subset ( (tTl-ti) *plot .unit , ti<=tTl & tTK=tf) # time stamps and loss rates 
y <- subset ( Tl, ti<=tTl & tTK=tf) 

ii <- round(seq(l, length(x) , length=2048) ) # choose a few representative points 
x <- x[ii]; y <- y[ii] 

plot(x, y, type="l", xlim=xl, axes=FALSE, xlab="", 

ylab="Loss rate [arb. units]") 
axis(2) ; box() 

x <- subset ( (tc-ti) *plot .unit , ti<=tc & tc<=tf) # collimator positions 

y <- subset(cp, ti<=tc & tc<=tf) 

plot(x, y, type="o", pch=20, xlim=xl, xlab="", 

ylab=expression(paste("Collimator position [" , mu*m, "]"))) 
mtextC'Time [min] " , side=l, line=2, outer=TRUE) 

#line 365 "dmcs.Rnw" 
# sample parameters 

Jci <- 0.040; Jcf <- 0.050 # initial and final collimator actions [urn] 

dt <- 1 # time to move collimator from Jci to Jcf [s] 

Ai <- 1 # initial slope of distribution function [um~-2] 

Af <- Ai*(Jci/Jcf) # final slope 

tt <- seq(0, dt, length=5) # choose times 

JJ <- Jc(tt, Jci=Jci, Jcf=Jcf, dt=dt) # calculate coll. pos. 
AA <- Ac(tt, Ai=Ai, Af=Af, dt=dt) # calculate slopes 

co <- gray (seq(0 . 6 , 0, length=length(tt) ) ) # define colors, grayscale 
Jmin <- 0; Jmax <- max(Jci,Jcf) # plot range 
plot ( c ( Jmin , Jmax) , 

c(0, max(Jci*Ai, Jcf*Af )) , 

type="n" , 

xlab=expression(paste("Action [" , mu, "m] " , sep="")), 

ylab=expression(paste("Distribution function [" , mu*m~{-l}, "] " , sep=""))) 
abline(v=JJ, col=co, lwd=2) # plot collimator positions 
for(i in seq(along=tt) ) { # plot asymtotic functions 

curve(flin(x, A=AA[i], J0=JJ [i] ) , Jmin, Jmax, add=TRUE, col=co[i], lwd=2) } 
legend ("bottomlef t" , bty="n" , # describe times for each curve 
title="Time [s] " , title . col="black" , 
legend=f ormat (tt , di=2) , col=co, text.col=co) 
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xx <- . 5* ( Jmin+Jmax) # labels curves 

text(xx, flin(xx, A=Ai, JO=Jci) , expression(f [0] ( J) ) , pos=2) 
text(xx, flin(xx, A=Af , J0=Jcf ) , expression(f [infinity] (J) ) , pos=4) 
xx <- 0.6*Jmax; yy <- . l*max( Jci*Ai , Jcf*Af) 
text(xx, yy, expression(f [a] ( J ,t) ) , pos=2) 
xx2 <- 0.7* Jmax 

lines (x=c(xx, xx2) , y=c(yy, flin(xx2, A=AA[2], J0=JJ[2]))) 
caption. text <- "no" # text for LaTeX caption 
if(Jcf<Jci) caption. text <- "inward" 
if(Jcf>Jci) caption. text <- "outward" 

#line 554 "dmcs.Rnw" 

# sample parameters 

Jl <- 0.040; J2 <- 0.050 # initial and final collimator actions [urn] 

dt <- 1 # time to move collimator from Jci to Jcf [s] 

Al <- 1 # initial slope of distribution function [unT-2] 

A2 <- A1*(J1/J2) # final slope 

D <- le-5 # diffusion coefficient [um~2/s] 

# choose times to plot 
tscale. short <- abs ( J1-J2) ~2/D 
tscale.long <- min( Jl , J2) ~2/D 

tt <- c(seq(0, dt , length=5) , # collimator moving 

c(3*dt, tscale . short , 2*tscale . short) , # short term after 
c(5*tscale . short , tscale.long, 2*tscale . long) ) # long term after 
co <- gray(seq(0.6, 0, length=length(tt) ) ) # define colors, grayscale 
Jmin <- 0; Jmax <- max(Jl,J2) # plot range 
layout (matrix (1 : 2, nrow=l, ncol=2)) # plot layout 
for(i in 1:2) { 

if(i==l){ # left plot, inward 

Jci<-J2; Jcf<-Jl; Ai<-A2; Af<-Al> 
if(i==2){ # right plot, outward 

Jci<-Jl; Jcf<-J2; AK-A1 ; Af<-A2} 
JJ <- Jc(tt, Jci=Jci, Jcf=Jcf, dt=dt) # calculate coll. pos . 
AA <- Ac(tt, Ai=Ai, Af=Af, dt=dt) # calculate slopes 

ylabel <- expression(paste( "Distribution function [" , mu*m~-[-l}, "]", sep="")) 
if(i==2) ylabel <- "" 
plot ( c ( Jmin , Jmax) , 

c(0, max(Jci*Ai, Jcf *Af )) , 

type="n" , 



xlab=expression(paste("Action [", mu, "m] " , sep="")), 
ylab=ylabel) 

abline(v=JJ, col=co, lwd=2) # plot collimator positions 
for(i in seq(along=tt) ) { # plot asymptotic functions 

curve(fv(J=x, t=tt [i] , Ai=Ai, Af=Af, Jci=Jci, Jcf=Jcf, D=D, dt=dt) , 
Jmin, Jmax, add=TRUE, col=co[i], lwd=2) } 
legendO'bottomlef t" , bty="n" , # describe times for each curve 
title="Time [s] " , title . col="black" , 

legend=f ormatC(tt , dropO=TRUE) , col=co, text.col=co, cex=0.8) 

} 

#line 668 "dmcs.Rnw" 
# sample parameters 

Jl <- 0.040; 32 <- 0.050 # initial and final collimator actions [urn] 
dt <- 1 # time to move collimator from Jci to Jcf [s] 

D <- c(le-6, le-5, le-4) # show the effect of different diffusion coeff 
k <- 1; B <- 

tscale. short <- abs ( J1-J2) ~2/min(D) # estimate time scales 
tscale.long <- min( Jl , J2) ~2/min(D) 
tmax <- tscale . short/2 

co <- gray(seq(0.5, 0, length=length(D) ) ) # define colors, grayscale 
Jmin <- 0; Jmax <- max(Jl,J2) # plot range 
layout (matrix (1 : 2, nrow=l, ncol=2)) # plot layout 
for(i in 1:2) { 

if(i==l){ # left plot, inward 
Jci<-J2; Jcf<-Jl 

ylabel <- "Loss rate [arb. units]" 
legend. pos <- "topright"} 
if(i==2){ # right plot, outward 
Jci<-Jl; Jcf<-J2 
ylabel <- "" 

legend. pos <- "bottomright"} 
ADD <- FALSE 
for(j in 1: length (D)){ 

Ai <- 1/D[j] # scale slope to keep steady-state rates constant 

Af <- Ai*Jci/Jcf 

if(j>l) ADD <- TRUE 

curve (loss . rate (t=x, Jci=Jci, Jcf=Jcf, Ai=Ai, Af=Af, D=D[j], 
dt=dt, k=k, B=B) , -0.1*tmax, tmax, 
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xlab="Time [s] " , ylab=ylabel, add=ADD, lwd=2, col=co[j])} 
if(i==l) legend (legend. pos, bty="n" , # print diffusion rate for each curve 
title=expression(paste("Dif fusion rate [" , mu*m~2/s, "]")), 
title. col="black", 

legend=f ormatC(D, format="E", di=l) , col=co, text.col=co) 
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